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Diffusion  First  Passage  Times:  Approximations  and  Related  Differential  Equations 


By 

Michael  L.  Wenocur 
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Palo  Alto,  California 


A  finite  spectral  expansion  is  presented  for  the  distribution  of  first  passage  to  a  fixed 
level,  by  a  diffusion  process  with  reflecting  lower  boundary.  The  n-term  expansion 
derived  here  matches  the  first  n  moments  of  the  passage  time  distribution.  We  derive 
also  an  interesting  representation  for  the  moment  generating  function  of  the  first 
passage  distribution. 
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1.  Introduction  and  Motivation 


Let  {ft  ,  F  ,  P)  be  a  probability  space.  Suppose  that  Ft  is  a  filtration  of  F 
and  that  (W(£),  t  >  0}  is  standard  Brownian  motion  adapted  to  Ft. 

Let 

A  =  +  n{x)D  (1.1) 

be  the  infinitesimal  generator  of  a  diffusion  X(t)  on  [0  ,  r]  satisfying 
dX(t)  -  a{x)dW{t)  +  n(x)dt  , 

with  a  reflecting  boundary  at  0  and  absorption  at  r  <  oo  ,  and  where  cr2( x )  >  0 
and  n(x)  are  continuously  differentiable  on  [0  ,  r]  . 

Define  the  stopping  time  rr  by 

Tr  =  inf  (s:  «>0  ,  X(s)  >  r} 

and  the  moment  generating  function  !f>(a :,y)  by 

*#(*,!/)  -  £{e'r’  |JT(0)-*}  s  £*(e'r’l* 

First  Passage  Times  as  Failure  Times 

Our  motivation  for  studying  first  passage  time  distributions  is  their 
relevance  to  modeling  of  failure  times.  Indeed,  this  paper  continues  the  line  of 
development  initiated  in  [9],  where  a  stochastic  process  is  used  to  model  system 
state,  ie,  wear- and- tear,  and  failure  occurs  when  either  a  traumatic  killing  event 
occurs  (killing  events  happen  with  rate  k(x)  in  state  x),  or  the  system  is  retired 
when  i vear-and-tear  reaches  some  predefined  threshold  (ie,  a  first  passage 
occurs). 

For  example,  if  system  state  is  modeled  as  Brownian  motion  with  positive 
drift,  then  first  passage  to  a  specified  threshhold  has  an  inverse  Gaussian 
distribution.  This  first  passage  distribution  has  been  successfully  applied  to 
numerous  problems  to  obtain  good  fits,  (cf  Jorgensen  [5]  ). 

A  related  but  parallel  line  of  development  is  explored  in  Wenocur  [12]  , 
where  the  killing  time  distribution  of  Brownian  motion  with  quadratic 
extinction  rate  is  calculated. 

Our  aim  in  this  paper  is  to  study  first  passage  time  distributions,  where  the 
system  state  process  is  a  general  diffusion  with  reflection  at  the  origin  and 
absorption  at  r  <  oo.  That  is,  the  system  state  evolves  as  a  diffusion,  and 
failure  occurs  at  the  epoch  of  first  passage  (or  absorption)  to  level  r  . 


In  future  work,  we  intend  to  explore  the  practical  ramifications  of  employing 
the  computational  methods  suggested  here  to  evaluate  interesting  first  passage 
times  statistics. 


Backward  Equation  for  First  Passage  Time  Distribution 

Let  w(x,t)  denote  the  tail  of  the  first  passage  time  distribution,  ie, 

w(x,t)  =  Pz{rr  >  t}  . 


The  backward  differential  equation  for  w(x,t)  is 


dw(x,t) 

dt 


M(x) 


dw(x,t) 

dx 


o^(i)  d2t v(x,t) 
2  dx2 


Aw(x,t) 


(1.2) 


for  (x,<)  e  (0,r)  X  (0,oo)  ,  with  boundary  conditions 

tt>(x,0)  ■■  1  for  0  <  x  <  r  ,  and  for  all  t  >  0  w(r,t)  ™0  and  »  o  . 

For  a  derivation  of  this  equation  and  other  related  quantities  see  [6,  pp  222- 
224]. 

The  Spectral  Representation  for  w(x,t) 

The  following  representation  for  w(x,t)  is  valid  whenever  cr^x)  and  p(x)  are 
sufficiently  smooth: 

M*. 0-  S  ck  l1-3) 

*-l 


where  ak  and  <f>k  are  eigenvalues  and  eigenfunctions,  and  ck  are  generalized 
Fourier  coefficients,  all  defined  below  (This  representation  is  proved  in  Section 
8). 

The  {  4>k  ,  k  >  1}  are  eigenfunctions  of  A  corresponding  to  the  eigenvalues 

{Qk  ,  *  >  l},  ie* 

A<t>k  —  —ak  (j)k  , 
and 


c*  -  /  <t>k{x)p{x)dx  , 
o 


where  p(x)  is  given  by 


and  7r(x)  is  given  by 


I 

7r(x)  =*  exp J  2p(u)/cr2(u)du  . 
o 


(1.5) 


In  general  an  arbitrary  function  /  6  L2(p)  will  have  a  Fourier  type 
expansion,  ie, 

/  =  £ck(f>k 

k- 1 

where  equality  is  interpreted  in  the  L2(p)  sense  and 

r 

Ck  -  I f  (x)<t>k(x)p(x)dx  . 

0 


Remark:  In  the  sequel  it  is  assumed  that  A ’s  eigenvalues  form  a  complete  set  in 
L2(p).  The  completeness  of  A' s  eigenfunctions  can  be  assured  by  certain 
regularity  conditions  on  the  infinitesimal  parameters  cr2( x)  and  p(x)  .  For 
example,  o2(x)>0  and  the  continuity  of  o2  (x)  and  p  (x)  are  sufficient 
conditions.  See  [11,  chap  lj  for  more  details. 


A  Generalization 


lV  ‘o  '* 


This  paper  is  primarily  concerned'  with  computing  first  passage  time 
statistics.  In  alluded/ to  in  (1.1);  a"  general  reliability  model  was  proposed 

in  which  system  failures  o*cur  when  either  system  wear-and-tear  reaches  some 
maximum  permissible  level  (ie,  a  first  passage  occurs)  ,  or  when  some  killing 
event  happens  (such  killing  events  occur  withrate  k(x)  in  state  x).  Under  this 
model  w(x,t)  satisfies  the-foHowing  equation. 

-  t, +  ?t»)  -  Bw[x,t)  , 


with  the  same  boundary  conditions  as  (1.2)  ,  and  where 

Bf  (*)  =  jo 2(x)}  "(x)  +  M*)/  (*)  +  *(*)/(*)  ■ 

>It  is  possible  to  solve  for  w(x,t)  and  related  quantities  with  methods  very 
similar  to  those  presented  here.  , 

s.. 

In  Section  2,  algorithms  for  approximating  w(x,t)  are  obtained.  In 
particular,  the  infinite  spectral  expansion  for  w(x,t)  is  approximated  by  an  n- 
term  sub-expansion  which  matches  the  first  n—  1  moments.  Section  2  concludes 
with  some  remarks  about  our  preliminary  computational  experience. 


Proofs  validating  the  spectral  expansion  and  the  related  approximation 
scheme  are  given  in  the  Appendix,  Section  8  of  this  paper. 

**  In  Sections  3  and  4,  methods  are  given  for  obtaining  the  eigenvalues  and 
first  passage  moments,  necessary  for  computing  approximations  to  w{x,t )  .  In 
Section  5,  computational  issues  related  to  calculating  the  moment  generating 
function  are  considered. 

S'^tions  6  and  7  include  theoretical  complements  about  first  passage  times. 
In  particular,  the  moment  generating  function  is  shown  to  possess  an  interesting 
representation  having  exponential  form  Jfef- -equation  (7.1)  )/This  exponential 
representation  is  related  to  asymptotic  expansions  used  in  analyzing 
perturbations  of  certain  second-order  differential  equations,  m  _ 
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2.  Approximating  The  First  Passage  Time  Distribution 

The  preceding  discussion  might  suggest  that  solving  for  w(x,t)  is  fairly 
straightforward.  Generally,  eigenvalues  and  eigenfunctions  are  difficult  to 
obtain.  However  the  problem  of  approximating  w(x,t )  can  be  approached  by 
the  method  of  moments.  One  technique  is  to  calculate  the  first  three  moments, 
and  then  use  the  Pearson  curve  fitting  method  (cf  [lOj).  This  method  is 
computationally  feasible,  and  the  Pearson  family  of  curves  includes  some 
important  first  passage  distributions,  such  as  the  gamma  distribution  (cf  [l]). 
The  merits  of  this  approach  will  be  studied  in  a  forthcoming  paper. 

Given  n  Eigenvalues  And  n  —1  First  Passage  Moments 

A  more  computationally  intensive  approach,  but  one  founded  on  stronger 
theoretical  grounds,  is  the  following.  Suppose  that  n  moments 
{Mk(x,r)  ,  1  <  k  <  n  }  are  known,  where  Mk(x,r)  —  Ex[  r*  ],  as  well  as  the 
first  n  eigenvalues  { ak  ,  1  <  k  <  n  }  .  Then  use  a  finite  sum  in  place  of  the 
infinite  sum  in  equation  (1.3).  In  particular,  approximate  w(x,t)  by 

wn(xit)  ~  E  P]'Xx) 

i- 1 

where  p^",  s  {p^  ,  ■  ■  ■  ,  p„w)  satisfies  for  0  <  k  <  n—  l 

«  .  Mk(x,r) 

E  pj  (*)  M  - - -r. -  where  Hj  =»  1  /a,  .  (2.1) 

i- 1  *• 

Ideally  we  want  wn(x,t)  to  be  a  distribution  function,  ie,  tun(x,t)>0  and 
iun(x,s+t)  >  wn(x,t)  whenever  s  >  0.  It  is  not  clear  that  solving  (2.1)  always 
produces  such  a  function.  This  issue  requires  further  investigation. 

Obtaining  The  Weighting  Factors 

The  weighting  factors  {pkn)  ,  1  <  k  <  n  }  in  (2.1)  above  are  obtained  by 
solving  the  following  linear  system: 


where  mk  «-  Mk{x,r)/k\  . 

Observe  that  the  matrix  {/** ,  1  <  j  <  n  ,  0  <  k  <n}is  none  other  than  the 


transpose  of  the  celebrated  Vandermonde  matrix.  Cramer’s  equation  gives  the 
following  formula  for  . 


Pkn)  =  E  '  •  •  ,P*-i  ,  P*+i,  '  '•  P»)  /  fl  {Hj  ~  Hk)  (2*3) 

j  -  1  j  —  I 

]+k 

where  gT  are  the  signed  symmetric  functions  defined  as  follows: 

<7o(al’fl2»  '  '  »am)  =  1 

and  for  r  >1 

?r(a1,a2,  •  •  •  ,am)  =  •  '  ’  «i,(-l)r  • 

1  <  i'i  <»'*<••  •  < <  m 

In  Section  8  the  following  convergence  theorem  is  proved: 

Theorem 

For  fixed  k 

I  P*(,)  -  ck<}>k{x) I  "  o(~V)  i 

n 4 

where  ck<f>k(x)  is  defined  by  equation  (1.3). 

Given  2n— 1  Moments  Only 

Suppose  that  the  first  2n—  1  moments  have  been  determined.  It  is  possible  to 
approximately  determine  the  first  n  eigenvalues  by  solving  the  following  system 


of  equations  for  /^,  ,  p/* 

,  i  < «  < 

n. 

■ 

1 

1 

1 

i 

i 

Mi 

M2 

Mn 

Mi2 

m22 

Mn 

- 

m2 

M2”-1 

m!*-1  •  ■  • 

M2”-1 

p.w 

where  >  0  . 


One  approach  is  to  solve  for  {  pjn )  ,  j=l,  ■  ,n  }  in  terms  of  {/ik  ,  1  <  k  <  n} 
and  {  m0  ,  m1  ,  •  •  •  ,mn_ j  }  as  detailed  above,  and  then  use  the  remaining  n 
constraints  to  determine  the  {  pk  ,  1  <  k  <  n  }.  This  reduced  system  can  then  be 


solved  using  mathematical  programming  techniques,  eg,  approximate  Newton- 
Raphson  techniques.  The  numerical  stability  and  feasibility  of  this  method 
merit  further  study. 


Preliminary  Computational  Experience 

Numerical  experiments  conducted  in  the  C  programming  language  suggest 
the  following.  For  small  values  of  t  ,  finite  difference  methods  are  more 
attractive  than  finite  spectral  expansions,  but  the  finite  expansion  approach  is 
preferable  for  large  values  of  t  .  For  relatively  small  number  of  terms  (less  than 
6),  the  method  of  moments  seems  to  be  an  attractive  method  of  computing 
weights,  especially  when  moments  are  also  to  be  computed.  Because  the 
Vandermonde  system  of  equations  grows  increasing  unstable  with  its  dimension 
(but  see  Bjork  [2]),  longer  expansions  require  that  the  spectral  coefficients 
should  be  computed  in  terms  of  the  eigenvectors  and  Fourier  coefficients. 

Small  values  of  t  will  generally  be  of  interest  for  first  passage  distributions 
with  small  means,  and  conversely.  Roughly  speaking,  if  /i(x)  >  0  or  a{x)  is 
large,  then  first  passages  occur  rapidly.  But  if  n(x)  «0  ,  or  n(x)  <  0  and  o(x) 
is  small,  then  first  passages  occur  slowly.  These  observations  should  help  to 
guide  the  choice  of  numerical  method. 

The  following  example  illustrates  the  efficacy  of  the  method  of  moments  for 
computing  the  spectral  weights.  Let  o(x )  s  1  ,  fi(x)  s  0  and  r  =  0,  ie,  driftless 
Brownian  motion  with  a  reflecting  boundary  at  0  and  an  absorbing  boundary  at 
1.  In  the  following  table  exact  and  approximate  values  for  ck<f)k( 0)  ,  k  —  1,...,  5 
(cf  eqn  (1.3)  with  x  —  0). 


Method  of  Moments 

Coefficient 

3  Terms 

4  Terms 

5  Terms 

Exact 

c^i(O) 

1.2731 

1.2732 

1.2732 

1.27324 

c  2^2(0) 

-0.4002 

-0.4226 

-0.4243 

-0.42441 

c  3^3(0) 

0.1272 

0.2228 

0.2500 

0.25464 

«  4^(0) 

-0.0735 

-0.1452 

-0.18189 

c  5^5(0) 

0.0463 

0.14147 

-7- 


3.  Solving  For  The  Eigenvalues 


The  Eigenvalue  Equation 

Before  the  eigenvalue  equation  can  be  introduced,  the  eigenfunction 


differential  equation  must  be  rewritten  in  more  suitable  form.  To  do  so,  let 
Yf(x)  satisfy 


^-o2(x)Y9"  (x)  +  n(x)Y9' (x)  +  dYg(x)  =  0 

4M 


(3.1) 


Multiplying  (3.1)  by  p{x)  gives 

7r {x)Ye{x)  +  ir'(x)Yf  (x)  +  9p(x)Y$(x)  , 

where  p(x)  and  tc{x)  are  given  by  (1.4)  and  (1.5)  respectively. 

Now  suppose  that  Yt  satisfies  the  boundary  conditions 
Yj( 0)  -  0  and  Y9( 0)  =  1  . 

Then  define  u^O)  by 

<40)  =  Y, (r)  . 


(3.2) 


(3.3) 


The  eigenvalues  of  equation  (3.1)  are  none  other  than  the  zeroes  of  u>  ,  see  [3, 
Chap  8]  for  further  details. 


Standardized  Eigenvalue  Problem 

It  is  also  possible  to  transform  (3.1)  into  more  standard  form  using  the 
following  transformation  scheme.  Setting 


du 


“((*.)/*)■/» 


reduces  (3.1)  to 


+  «*)4r  +  0y=o, 


dz‘ 


dz 


(3.4) 


where  (3{z)  =  {//(*(*))  -  ~ (*(*))}/(< t2(z(^))/2)1/2  . 

4 


Putting  Y(z)  =  g[z)y(z )  ,  where  g(z)  =  e' 

£jl 


gives 


dz 


2  +  (0  “  ?(*))?  =  0 


(3.5) 


-8- 


I 


where  q(z)  =  ■^■/32(z)  +  and  with  boundary  conditions  y  (0)  =  0  and 

y(b)  =  0  with  6  =  Z(r )  . 

The  eigenvalues  of  (3.5)  are  the  same  eigenvalues  as  those  of  (3.4)  and  (3.2)  . 
Moreover  the  eigenfunctions  of  (3.5)  are  easily  transformed  into  those  of  (3.2)  . 
In  particular  if  ipn{x)  1S  the  eigenvalue  corresponding  to  ocn  for  (3.5)  ,  then 

<t>n(x)  =  9(Z(x)Wn(Z(x)). 


The  following  asymptotic  results  (as  n  — ►  oo  )  are  known  about  the  eigenvalues 
and  eigenfunctions  of  the  standardized  problem  (cf  [11,  p.  19]): 

cxn  -  (n+l/2)V/ft2  +  0(1)  ;  (3.6) 

V^n(x)  =  (2/6)1^2cos((n+l/2)7rx/6)  +  0(— )  ;  and  (3.7) 

n 

V’n(i)  =  —  (n+l/2)7r(2/63)1/2sin(n7r;r/6)  +  0(1)  . 


i 

i 

v 

* 

4 


► 

) 

i 


(3.8) 


4.  Obtaining  First  Passage  Moments 

To  solve  (2.1)  we  need  to  produce  the  moment-sequence 
{Mjc(x,r)  ,  1  <  k  <  n},  three  such  methods  are  outlined  below. 


Complex  Integration  To  Invert  Moment  Generating  Function 

Corollary  6.4  shows  that  ^6(x,r)  is  an  analytic  function  of  9  around  0,  and 
so  Cauchy’s  formula  implies  the  identity 


Mn{x,r)  ^  j_  ^9{x,r) 

«!  27r*  \i\-e0  0n  +  l 

where  60  is  sufficiently  small. 


(4.1) 


The  integrals  in  (4.1)  may  be  computed  numerically  using  Gaussian 
quadrature  to  minimize  the  number  of  values  of  9  to  be  evaluated,  and  then 
taking  the  real  part.  This  reduces  evaluating  the  equation  to  calculating  a 
small  number  of  values  of  'frg(x,r)  .  Evaluating  &g(x,r)  may  be  done  by  either 
using  finite  differences  to  solve  the  boundary  value  problem  (cf  equation  (7.6)), 
or  by  using  the  series  method  suggested  in  the  differentiation  approach,  or  some 
hybrid  of  series  and  finite  differences.  An  important  virtue  of  estimates  of 
Mn(x,r )  based  on  formula  (4.1)  is  that  the  accuracy  of  these  estimates  is 
independent  of  the  accuracy  of  the  n— 1  smaller  moments,  unlike  the  methods 
given  below. 


Recursive  Integration 

This  approach  iteratively  uses  (6.1)  to  compute  successive  moments.  This  is 
feasible  when  the  successive  moments  form  a  closed  family  of  integrals  (compare 
example  1),  or  when  only  a  few  moments  are  desired. 

Example  1 

Choosing  parameters  o2( x )  =  2(x  +  a0)  and  /i(x)  =  v  where  v  j=-  0 
gives  rise  to  the  iteration: 

r  w 

Mn{x)  =  f  - — -f — —  /  +  <r0 Y~ldudw 

x  (w  +  <Xq)  0 

An  easy  induction  will  show  that  for  u  not  integer  Mn_l(x )  satisfies  the 
expansion 

Mn_l(x)  =  c[n, 0]  +  £  c[n,£](x+a)*  +  £  d\n,k\{x+a,)k~i' 

k- i  k-i 


where  the  coefficients  are  calculated  using  the  iteration: 


c[n,fc]  =  —nc[n—l,k—l]/(k(k—l)+vk)  for  k  >  1  ,  and 

d[n,k\  =  —nd[n—l,k—l\/((—i/+k)(—u+k—l)+v(—v+k))  for  k  >  2 

Finally  c[n,0]  and  d[n,l]  are  determined  by  solving  the  two-dimensional  linear 
system  arising  from  the  boundary  conditions  Mn(r)  =  0  and  Mn( 0)  =  0  . 

For  integer  u  ,  closed  formulas  for  all  moments  are  obtainable,  but  the 
calculations  will  be  messier. 

Differentiating  the  Moment  Generating  Function 

Successive  moments  may  be  obtained  by  calculating  the  ^-derivatives  of  the 
moment  generating  function  &g(x,r)  at  9  =  0.  This  approach  is  facilitated  by 
Kent's  observation  in  [8]  that  &g{x,r)  =  Tg(x)/  Tg(r)  where  Tg(x)  satisfies 

~<^(x)Tg  (*)  +  V{*)Tt{x)  +  0Tg{x)  =  0  (4.2) 

with  initial  conditions 

7,(0)  =  0  ,  T,( 0)*0.  (4.3) 

To  solve  for  the  ^-derivatives  of  #g(x,r)  it  suffices  to  solve  for  the 
derivatives  of  Tg(x).  We  can  obtain  Tg(x)  using  the  series  expansion  method 
around  r=0  to  solve  (4.2)  .  Under  certain  regularity  conditions,  the  Taylor 
series  coefficients  may  be  differentiated  with  respect  to  0  ,  and  the  series 
summed.  This  process  is  illustrated  in  Example  2  below. 

Example  2 

We  indicate  how  the  technique  in  (4.3)  may  be  applied  to  Example  1: 


Ux)  =  Z  bj{6)xi  . 
y-o 

Equation  (4.2)  implies  that 


<r(x)  Z 

j-o 


00  '>+  2 
2 


bj+2(0)x3  +  n{x)  Z  ti+l)b}+y(e)xJ  +  e  £  bj{9)x] 
y-o  y-o 


=  0 


Since  p(x)  =  u  and  a(j)  =  2(x-F<x0)  we  deduce  that 

9b j -2(0)  T  (^(j-l)U-2))bj-i(0)  +  <x0j(j~l)bj(9)  =  0  ,  for  j  >  2  , 


(4.4) 


rwv  wvjtj  wu  m  vuvu  hu  w  w.  i  rm  nwn*  rm  rjr  nx  *jl  rji  tut  hji  f j*  k  x*x  *.  if  wrKjni  URyRuunu/iun'iiT^^jiVTwyii^LNWiwm  wuwww  u 


6O(0)  =  1  and  b^O)  =  0  ,  for  j  =  0 , 1  . 


Repeatedly  differentiating  (4.4)  will  give  successive  iterative  formulas  for 
computing  the  Taylor  series  coefficients  of  T^n\x)  .  For  example,  if  n  =  1 


Obj-iiQ)  +  bj_2(6)  +  {v  +  (j — 2)(j — l))6yl_! [0)  +  Ob  y(y-l)&/(0)  =  0  ,  for  j  >  2  , 


with  initial  condition  given  by  equation  (4.4)  with  0  =  0. 

If  the  above  iteration  diverges,  we  can  always  try  renormalizing  by  x  and 
calculating  bj’,\0)xn .  If  r  is  sufficiently  small  then  renormalization  will  suffice, 
otherwise  Tg(x)  can  be  calculated  by  successively  moving  out  from  0  towards  r 
as  suggested  in  Section  5  below,  and  then  using  the  contour  integration  method 
given  in  the  beginning  of  this  section. 


5.  Some  Remarks  About  Computing  &9(x,r) 

Computing  ify(0 ,r) 

The  series  expansion  method  may  not  permit  solving  for  !fy(0 ,r)  in  a  single 
step.  However,  suppose  the  series  converges  for  some  y  6  (0,r)  ,  ie,  it  is  possible 
to  compute  &9{0,y)  by  the  series  method.  We  may  use  $9(0,y)  as  a  bootstrap 
to  calculate  !fy(y,r)  as  follows.  Observe  that  <fy(0,r)  =  ^(0 ,y)&${y,r)  .  Thus 


^${y,r)  -  - 


^(o.y)  d*t(y,r) 

d^fjO  ,y)  dy 

dy 


Using  this  initial  condition  and  Kent’s  normalization  technique,  it  is  possible  to 
calculate  V $(y  ,r )  starting  from  y  rather  than  from  0. 

Interpolating  ^(x,r)  And  A  Related  Boundary  Value  Problem 

Suppose  that  &9(x,r)  and  5fy(y,r)  have  been  obtained  ( x  <  y)  ,  and  it  is 
desired  to  calculate  V 9(z,r )  for  j  6  (i,y)  .  The  multiplicative  character  of 
&9(x,r)  implies  that 

It  thus  suffices  to  determine  V9{z,y)  .  We  have  ^9{x,y)  =  &9(x,r )/\fr9(y,r )  and 
—  1  .  Therefore  it  suffices  to  find  h9(z )  (  =  9(z,y ))  such  that 

■j-o2(z)h9" (z)  +  fi(z)h9  (z)  +  0h9(z )  =  0  (5.1) 

mi 

with  boundary  conditions  h9{x)  —  ty9(x,r)/ V9(y ,r )  and  h9(y)  —  1  .  These 
boundary  conditions  uniquely  determine  h9  .  To  solve  for  h9  ,  first  find  £0  and 
satisfying  (5.1)  ,  where  f,(a:)  =  1— t  and  ?,(x)  =*  t  ,  for  t  =  0,  1  .  Then  set 

h9{z)  -  ^9(x,y)^(z)  +  e0(^)(i  -  fy(z>y)£i(iO)/£o(sO  • 


6.  Theoretical  Complements 

In  this  section  some  of  the  properties  of  moments  of  tz  are  examined,  but 
first  some  new  notation  is  introduced. 


Define  Mn(x,y)  =  Ez{Ty)  for  x  <  y  and  n  >  0  ,  and  let 
dMn(x,y) 


Mn(x,y )  = 
Mn\x,y)  = 


dx 

d2Mn(x,y) 

dx2 


Recursive  Equations  For  Moments  Of  rx 

The  functions  Mn(x,y)  jointly  satisfy  the  iterative  differential  equation  (cf 
[6],  p.  203,  equation  (3.38)  )  : 

~<^(x)Mn(x,y)  +  fi(x)Mj{x,y)  +  nMn_x(x,y)  =  0  (6.1) 

subject  to  Mn(0,y)  =  0  and  Mn(y,y)  =  0  . 

Lipschitz  Conditions  For  Moments  Of  Tx 
Lemma 

Mn(x,y)  is  a  smooth  function  in  x  and  y  jointly  ,  and  there  exists  a  constant 
C  such  that 


\Mn(x,y)\<Cny2n  \y-x)n\  (6.2) 

and 

KW)l<CV(n_I)n! 

for  all  x  <  y  . 

Proof: 

Set 


«(*)  -  exp{-/-^^-  d£  } 


and 


m(x)  *=  \ /\a2(x)8(x)\  . 


Rewriting  (6.1)  as  follows 


d  Mn\x,y) 


dx  [  s(x ) 
implies  that 


=  —2nMn_1(x,y)m(x) 


Mn(x,y)  =  2 nj  [/M^^y )m(€)df]s(7?)dr? 


By  virtue  of  continuity  there  exists  a  constant  K  such  that 
II  *  II  s  nSup  |  s{x)  |  <  K 

0  <  x  <  r 


II  m  ||  <  K 


Therefore 


\Mn(x,y)\<2nj  {  K^WM^WiZdv 

x  0 


V 

jS 

'll 

*-r 


$ 


2n  A:2||Mn_1||j/(y-i)  <  2iC2y2n  || 


An  easy  induction  implies  that  Mn(x,y)  is  a  smooth  function  in  x  and  y  , 
moreover 


||A/„(z,y)||  <  2nA2ny2n-1(y  — x)n! 


||A/n(x,3/)||  <  2 nK2ny2^~%\ 

Taking  C  =  (2 K2)  completes  the  proof. 
(6.4)  Corollary 

^#(x,y)  <  oo  whenever  \#\  <  C  ly  2,  and 


*#(*.y)  -  E  Mnix,y) 


Infinitesimal  Relations  Governing  First  Passage  Moments 


E 


Proposition 

Define 


dMn(x,  y) 


dy 

and 

un  ix)  = 

Then  Mn{x,y)  ,  Un(x,y)  ,  un(x)  satisfy 


A/„(0  ,y)  =  £ 

j-  o 


/  \ 
n 

j 


Mj{0,x)Mn_j(x,y ) 


(6.5) 


^„(0,x)  =  V 

y-o 


/  \ 
n 


J 

V  / 


(6.6) 


Proof 

Conditional  on  X(0)  =  0  the  strong  Markov  property  (SMP)  implies  that 
{  ra,  »  ra2  —  ra,  *  •■•*Tr  —  ra,  }>  where  0  <  a,  <  o2  <  •  ■  •  <  o,  <r,  form  a  set  of 
independent  random  variables.  In  particular 


M,(0,y)  =  E0[rny\  =  £°((r,  +  ry  -  r,)"] 


n 

=  r 

y-o 


J 

V  / 


-  E 

y-o 


£°K)J  j  £*((rf-rx)->  ] 


which  coincides  with  (6.5)  .  Using  a  little  algebraic  manipulation  on  (6.5)  shows 
that 


n  — 1 


\Mn{0,y)-Mn{0,x)}/{y-x)  *=  £ 

y-o 


y  jMj(0,j)Mn_;(x,y)/(y-x) 
Now  letting  y—+x  in  (6.7)  yields  (6.6)  .  QED 


1.6-7) 
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Comments 

Equations  (6.5)  and  (6.6)  provide  some  nice  intuition  about  the  way  that 
first  passage  times  from  x  to  y  depend  on  first  passage  times  from  0  to  x. 

Equations  (6.5)  and  (6.6)  are  similar  to  (6.1)  ,  but  may  capture  better  the 
dependence  of  higher  moments  on  lower  moments.  From  a  practical  point  of 
view  equation  (6.1)  is  certainly  preferable  for  moment  calculation.  In  Section  4, 
other  methods  are  proposed  for  calculating  the  moments  Mn(x,y)  . 


7.  A  Representation  Result 


Theorem 


oo  y  Qn 

#t{x,y)  =  exp{  £  J  un(z )  —  dz} 

n  —  1  i 


11 ! 


where  {un(z)  ,  n>  1}  satisfy 

+  M*)«i(*)  =  1  < 
and  for  n  >  2 


^{x)un{x)  +  n(x)un(x)  -  jcr(i)  £ 


“*(*K-*(*) 


Proof 

We  begin  by  showing  that  (7.1)  holds  for  #f(Qj)  . 

Let  xy  »»  jr  /L  for  0  <  j  <  L .  Using  the  SMP  as  in  the  proof  of  (6.5), 

L-i 

W.O  *  77  1)  • 

j-0 

Taking  logarithms, 
i-i 

logtf,(0,r)«  27  log<^(x;,x>+1)  . 

>-  o 


Applying  the  proposition  in  Section  8  to  the  above  yields 

log ^#(0,r )  -  27  [#*{xj’Xj+ j)  -  1  +  {*t{Xj,xj+l)-l)2g(*,{xrxJ+i))]  . 

j-o 


Since  *#(*.*)  -  1  it  follows 
d*t 

+  o,L- )  L- 

where  0  <  a  <  1  .  Also  since  'Pt(x,y)  is  a  smooth  function  in  x  and  y 
d*pJz,y) 

it  follows  that  — -  ’  —  is  uniformly  bounded  on  the  region  0  <  x  < 
dy 

there  exists  a  function  h(L)  such  that  A(L)~~0(1) 
l°g^(0,r)  =  27  l*t(xj’xj+ i)_1l  +  h(L)L~l  . 


(7.1) 


(7.2) 


(7.3) 


jointly  , 
<  r.  So 


Replacing  each  lfy(xy,x;+1)  in  the  above  by  the  expression  in  (6.4),  yields 


log^(0,r) 


L- 1 


E 

3-0 


E 


n  — I 


OnMn{Xj,X y+i) 

n! 


+  0{L~'). 


Since  the  summands  are  positive  the  order  of  summation  may  be  permuted  to 
get 


oo  L- 1  dnMn(x:,Xj  . .)  , 

log#,(0,r)-  S  £  - +  0(L-'). 

n  —  1  /-  0  n- 

The  inner  summation  can  be  expressed  as  a  Riemann  sum 


log  ^,(0 


00  L- 1  Af.(x1,x,.1)  1  , 

»r )  -  E  &n  E  — \L; _r  -  T  +  )  • 


n_i  y- 0  n!  L‘ 

Equation  (6.2)  implies  that 

-  c”(x'+,)2*  -  c“ 
where  C0  *■  Cr2  . 


(7.4) 


The  Lebesgue  dominated  convergence  theorem  (applied  to  the  product  space 
{l,2,..}  X  (0, r ]  endowed  with  product  of  the  counting  measure  with  the 
Lebesgue  measure  on  [0,r])  implies  the  right  side  of  (7.4)  approaches  the  limit 

log*#(0,r)  -  £  ~-r  /  un{z)dz 

n-J  0 


as  L  — ►  00  ,  or 


^#(0,r)  -  exp{  27  “  /  «„(*)  *  }  • 

n-l  ”•  0 

Due  to  the  multiplicative  nature  of  V i(x,y )  it  is  easy  to  show  that 


W*,y)  -  exp{  £  /  *„(*)  —  dz  }  . 
«— 1  1 


(7.5) 


The  function  ^§{x,y)  satisfies  the  following  differential  equation  (cf  [6,  pp  203)  ) 

1  ,  d*#(x,y)  .  v 

-<r(i) - — - +  l(x) - — - +  00>(x,y)  -  0 


5x2 


ax 


(7.6) 


s 


dVf(0,y) 

subject  to  - - - =  0  and  'I'giy, y)  —  1  . 

ox 


i 


Substituting  equation  (7.5)  in  (7.6)  ,  shows  that  the  exponent  in  (7.5)  satisfies 
the  differential  equation 


1  „  oo  a*  „  oo  fl*  .  oo  an 

7»2(*)[(  E  t »«(*))  ~  E  ~7«n(^))  -  »{*)  E  ~Tun(*)]  +  0  =  o 

2  n-1  n!  n— 1  n!  n— 1  n! 


Since  the  above  series  converge  absolutely  for  #  sufficiently  small,  we  can 
rearrange  terms  to  obtain  a  single  power  series  in  9  .  Because  this  power  series 
is  zero  for  9  sufficiently  small,  all  its  coefficients  must  be  zero.  Setting  the 
coefficients  of  9  to  zero  yields  equations  (7.2)  and  (7.3)  . 

The  initial  condition  that  =  o  in  (7.6)  implies  un(0)  =  0,  where  n  >  1  . 

OX 

This  concludes  the  proof. 


s 
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s 
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£ 
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8.  Appendix 

Spectral  Representations  For  First  Passage  Time  Distributions 

Theorem  The  right-hand-side  of  (1.3)  is  the  unique  function  satisfying  (1.2), 
jointly  continuous  in  x  and  t  on  [0,rj  X  (0,oo)  with  absolutely 

integrable  over  [0,r]  X  ( N~l ,  N)  for  all  N  >  1. 


Proof: 

Suppose  that  w(x,t)  satisfies  equation  (1.2)  and  the  integrability  conditions. 
Observe  that  for  fixed  t  >0  w(x,t)  is  a  continuous  function  of  x  belonging  to 
L2(p)  (  where  p  is  defined  by  equation  (1.4)  ).  Therefore  w(x,t)  possesses  the 
orthogonal  expansion 


OO 


w(x,t)  -  £  ck{t)  <f>k(x) 

*-i 


where 


ck(t)  "*  /  w(a? ,t)  <t>k(x)p(x)dx  . 
o 


Multiply  (1.2)  by  <t>k{x)  p(x)  and  integrate  over  [0,r]  to  get 
/  dW <f>k(x)p(x)dx  >•  J  Aw(x,t)<j>k{x)p(x)dx  . 


I  ,  .  „ 


Now  since  Af  (x)p(x)  =  7r(x)/  (x)  +  tt  (x)f  (x)  (cf  (3.1)  and  (3.2)),  a  simple 
integration  by  parts  shows 


r  r 

f  Aw(x,t)<f>k(x)p(x)dx  =J  w(x,t)A  <j>k(x)p(x)dx  . 
0  0 


The  relationship  A<t>k{x)  —  —ak4>k(x)  implies  that 
/  4>k{x)p{x)dx  =  -  Jakw{x,t)<t>k{x)p{x)dx  . 


Integrating  both  sides  with  respect  to  t  over  [u0,u]  an<^  permuting  the  order  of 
integration  on  the  left-hand-side  yields  (permissible  because  Fubini’s  Theorem 

which  is  absolutely  integrable) 


/  /  -W ^ ? ^-1  <f>k (x )p(x )  dt  dx  =  /  /  —akw(x,t)<f>k(x)p(x)  dx  dt  . 
at 


o  «0 


o 
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Using  the  definition  of  ck(t )  on  the  the  last  equation  implies 


ck(u)  +  C  =  f  -Qkck(t)dt  , 


«0 


where  C  is  an  arbitrary  constant. 


Therefore  ck(t)  =  cke~aii .  It  remains  to  determine  the  constants  ck  ,  but  they 
may  be  derived  from  the  boundary  condition  w(x, 0)  =  1  for  0  <  x  <  r  as 
follows.  Since  p  is  continuous,  u>(:r,0)  =  1  for  almost  all  p{dx),  and  1  e  L2(p). 
So 


oo 

E  ckMx) 

k- 1 


(8.1) 


where 


ck  =  J  <f>k{x)P(x)dx  . 
o 


(8.2) 


Now  Theorem  1.9  of  [11]  implies  that  the  right-hand-side  of  (8.1)  converges 
pointwise  to  1  on  (0 ,r)  .  Therefore  w(x,t)  has  the  representation 


w(x,t)  =  £  ck  e  ‘  <j>k{x)  . 


(8.3) 


To  prove  the  converse,  suppose  that  w(x,t)  is  defined  by  (8.2)  and  (8.3)  jointly. 
Equations  (3.7)  and  (8.2)  imply  that  the  coefficients  ck  ,  k>  1  are  uniformly 
bounded.  Hence  for  t  >  e  >  0  the  series  converges  uniformly  to  a  function 
continuous  on  the  product  [0,r]  X  [e,oo)  .  The  uniform  convergence  and 
boundary  conditions  on  the  eigenfunctions  imply  the  boundary  conditions  on  x 

.  The  integrability  conditions  on  ^(-uO  f0||0W  in  similar  fashion.  The 


dt 


boundary  condition  w(x,0)  =  1  for  x  £  (0,r)  follows  from  Theorem  1.9  of  [11] 
and  equation  (8.1)  The  differential  equation  (1.2)  may  be  derived  from  the 
definition  of  derivatives  as  limits  of  divided  differences,  and  the  dominated 
convergence  theorem  applied  to  series.  QED. 


It  should  be  noted  that  the  first  passage  time  distribution  satisfies  the 
regularity  conditions  of  the  theorem,  and  therefore  must  have  representation 
(1.2). 


Convergence  Of  The  Finite  Approximations  To  The  Infinite  Vector 


It  will  now  be  shown  that  the  solution  vector  to  system  (2.1)  ,  denoted  by 
p(,)  ,  l  <  j  <  «},  converges  at  rate  n-2  component-wise  to  the  infinite 

vector  p  =  { pj  ,  j  >1}  ,  where  =  Cj<f>j(x)  . 


Theorem 

I  p*1;1  -  p*0  I  =  0(-r) 

n 

Proof:  Define  ul']  ,  S(k']  and  //**’  as  follows: 

n  , 

i^i*1  =  27  t^jPj  where  0  <  £  <  n— 1 
j-i 

=  mk  —  i/**’  =  27  VjPj  where  0  <  k  <  n— 1  (8.4) 

j-n  +  i 

Vk  ~Pk  Pk 

Observe  that  ,  1  <  A:  <  n}  is  the  solution  to  (2.2)  where  {m*  ,  1  <  ^  <  n} 
has  been  replaced  by  {6^"’  ,  1  <  k  <  n}.  In  particular 

t?*"’  =  27  hi%)gn_j(n up2,  ■  ■  ■  ,Uk-\  .  ^*+1.  ■  P»)  /  ]1  (V]  —  Vk)  •  (8-5) 

1  <  y  <  n  1  <  y  <  n 

;>r 


For  1  <  7  <  n 

|  gn-j((p  1,^2-  •  •  -  ^*+i.  •  •  •  /O  )  I  <  Vk~j~k  • 

Also  equations  (3.2),  (3.6)  and  (3.7)  together  imply  that  |  c„  |  =  0(— )  ,  thus 

n 

(8.4)  and  (3.6)  imply 
|  6}*>  |  <  Cn"2'  . 

So 

I  27  ffn- M#2,  ■  ■  ■  #k-x .  Pk+u  -p.))  6}"  I  <  27  »rj~k I  Vn)  I 
y-i  y- i 

<  i;  =  rf-‘-'o(-r)  • 

i-i  n 


Thus 

I  £  9n-j((Pi’P 2-  ■  -  P*+I.  •  •  ■  0.)  )hjm)l  =  0( 

y-i 


n— *  — 1 


n 


Finally,  the  denominator  of  (8.5)  may  be  written  as 


n  (My  -  Mi)  =  Mi  *4(n) 


i  <  y  < « 

y»*c 


where 


4n)=  n  (My -Mi)  n 

i  <  y  <  *-i  *+i  <  j 

Now  equation  (3.6)  and  the  Weierstrass  Product  Convergence  test  jointly  imply 
that 

lim  din^  =  d  >  0  . 

n  —*>oo 

Thus  setting  k  =  kg  yields  rjj^  —  0(  £■)  as  claimed. 


A  Conjecture 

It  is  interesting  to  note  that  pk(x)  is  linearly  proportional  to  the 
eigenfunction  <f>k(x),  and  therefore  Apk{x)  —  akpk{x).  This  suggests  that  pk'{x) 
will  approximately  satisfy  this  relation.  Observe  that  Amk( x)  = 

Thus 

Apk*\x)  =  — my </n_y_i(^i,M2>  '  '  iPk-i  >  lik+ 1>  '  '  #*»)  /  11  (Mj  M*) 

i<y<n-i  1  n 

Comparing  this  with  (2.3)  suggests  that 

i((Mi,Ms,  •  •  •  ^*-1  .  ^*+1. '  '  ‘  /*»)  )  _  _  J_ 

lim  - 7" - ‘  TC'  ~  ak  — 

n— *oo  <7n_y((/* '  '  i^i-i  >  J  Mi 

Proposition 

For  1  <  x  <2 

|  /o?(x)  -  x  +  1 1  =  (x-\fg{x)/2  ,  where  |  <7(1)  |  <  1  • 

Proof:  The  mean  value  theorem  applied  to  logix )  at  x=l  gives: 
log(x)  —  log(l+x—l)  =Iog(l )  +  (x—  1)  —  (x— 1)  /2(l+cx(x— l)  ) 

where  0  <  a  <  1  .  The  prop  now  follows  from  x  >  1  and  c*  >  0  . 
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